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Abstract 

This  paper  reports  on  the  development  and  numerical  implementation  of  a  comprehensive  3D  computational  fuel  cell  dynamics  code  for 
PEMFCs.  The  code  solves  a  set  of  coupled  non-linear  conservation  equations  (mass,  momentum,  species,  energy,  electrical  potential  and  liquid 
water  saturation)  for  an  entire  unit  cell.  A  phenomenological  model  for  water  transport  in  the  membrane  is  solved  separately  for  the  membrane 
domain,  in  conjunction  with  calculation  of  the  water  content  on  the  boundary  such  that  that  water  balance  is  satisfied  on  both  sides  of  the  membrane 
interface,  and  the  numerical  implementation  of  the  model  is  validated  against  an  analytical  solution.  The  global  polarization  curve  predicted  with 
the  CFD  code  is  found  to  compare  favorably  with  reported  data.  A  detailed  validation  of  the  CFD  code  against  spatially  resolved  experimental  data 
is  presented  in  a  companion  Part  2  paper,  and  in  this  paper  base  case  simulations  for  a  unit  cell  with  straight  channels  are  presented  to  illustrate 
and  analyze  basic  physical  features,  transport  of  species  along  the  channel  and  coupling  between  heat  and  mass  transfer  processes.  Analysis  of 
the  results  shows  that  many  of  the  variables  of  interest,  including  mass  fractions  and  current  densities,  exhibit  similar  profiles  along  the  channel, 
which  suggests  that  reduced  dimensional  model  based  on  appropriate  similarity  variables  might  be  suitable  for  rapid  calculations. 

©  2008  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

The  last  few  years  have  seen  significant  developments 
in  computational  fuel  cell  engineering  (CFCE)  allowing 
multi-dimensional  simulations  of  coupled  transport  in  proton- 
exchange  membrane  fuel  cells  (PEMFCs),  e.g.  [1-9].  CFCE 
tools  can  provide  invaluable  assistance  in  analyzing  thermal 
and  water  management  problems  in  a  fuel  cell,  in  design  and 
optimization,  in  guiding  experimental  investigations,  and  ulti¬ 
mately  in  improving  performance  and  achieving  cost  reductions. 
But  many  challenges  remain  [10].  One  of  the  hurdles  in  firmly 
establishing  the  reliability  of  CFCE  tools  -  a  prerequisite  for 
their  systematic  use  in  real  design  -  is  the  lack  of  validation.  In 
principle,  a  numerical  solution  should  be  validated  against  exper¬ 
imental  data  of  equal  dimensions  for  all  the  variables  solved  in 
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the  computation.  This  is  impractical  for  most  fuel  cells  because 
of  the  limitations  in  existing  experimental  techniques  and  the 
inherent  difficulties  of  in  situ  measurements. 

State-of-the-art  multi-dimensional  CFCE  tools  are  typically 
built  on  the  well-established  computational  fluid  dynamics 
(CFD)  framework,  which  provides  capabilities  for  numerical 
discretization  and  solution  of  coupled,  non-linear  convection- 
diffusion  equations  that  govern  a  variety  of  thermo-fluid 
processes.  In  addition  to  the  conventional  CFD  framework, 
the  CFCE  tools  often  involve  considerations  of  transport  pro¬ 
cesses  specific  to  PEMFCs,  including  multi-phase  flows  in 
porous  media  and  microchannels,  electrochemical  reactions, 
and  transport  of  charged  species  and  water  in  ionomer  phase. 
Although  most  current  commercial  CFD  codes  can  provide 
numerically  robust  and  physically  reliable  simulations  of  com¬ 
mon  thermo-fluid  problems,  their  capabilities  in  handling  the 
additional  complexities  of  a  PEMFC  problem  still  require  great 
care  to  ensure  robustness.  Furthermore,  detailed  assessment  and 
scrutiny  of  some  of  the  physical  models  implemented  in  cur- 
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Nomenclature 

a  water  activity 

B  body  force  (N) 

Co2  oxygen  concentration  (molm-3) 

Cq2  oxygen  concentration  in  GDL  at  cathode  catalyst 
layer  (mol  m-3) 

Cq2  cathode  reference  water  concentration  (mol  m-3) 
Cp  specific  heat  ( J  mol- 1  K- 1 ) 

Cpi  specific  heat  for  ith  species  (J  mol-1  K-1) 

D[  diffusion  coefficient  for  ith  species  (m2  s-1) 

D'  diffusion  coefficient  (m2  s-1) 

Dx  water  diffusion  coefficient  (mol  m- 1  s  - 1 ) 

E°  equilibrium  cell  voltage 

F  Faraday  constant,  96487C 

g  Gibbs  free  energy  (J  kg- 1 ) 

gi  Gibbs  free  energy  for  species  i  (J  kg-1) 
h  mixture  enthalpy  ( J  kg  - 1 ) 

hi  enthalpy  of  ith  species  (Jkg-1) 

hi  sensible  heat  of  ith  species  (Jkg-1) 

i  current  density  (A  m-2) 

J  mass  flux  (kg  m-2  s-1) 

jo  exchange  current  density  (A  cm-2) 

M  molecular  weight  (kg  mol  - 1 ) 

Mm  equivalent  weight  of  a  dry  membrane  (kg  mol- 1 ) 
rhi  mass  flow  rate  of  species  i  (kg  s-1) 

rh\  phase  change  rate  (kg  s - 1 ) 

Nq  number  of  gas-phase  species 

rid  electro-osmotic  drag  coefficient 

p  pressure  (Pa) 

q  heat  flux  (W  m-2) 

R  membrane  resistance  (£?) 

Ru  universal  gas  constant  (8.314  J  mol-1  K-1) 

RH  relative  humidity 

5h  enthalpy  source  due  to  phase  change  (W  m3) 
(S/V)eff  effective  surface  to  volume  ratio  (m2  m3) 
s  saturation 

T  temperature  (K) 

t  time  (s) 

u  bulk  fluid  velocity  (ms-1) 

V  cell  voltage  (V) 

v  X-coordinate  (between  channel  and  land  area) 

y  7-coordinate  (perpendicular  to  MEA) 

Yi  mass  fraction  of  ith  species 

z  Z-coordinate  (axial) 

Greek  symbols 

A h®  enthalpy  of  formation  (J  mol-1) 

a  transfer  coefficient 

8  average  pore  size  (m) 

s  wet  porosity 

y)  activation  overpotential  (V) 

k  thermal  conductivity  (W  m-1  K-1) 

k  permeability  (m2) 

X  water  content 


/x  dynamic  viscosity  (kg  m-1  s-1) 
p  mass  density  of  mixture  (kg  m-3) 

Pdry  density  of  dry  membrane  (kg  m-3) 
pm  density  of  a  dry  membrane  (kg  m-3) 
a  electrical  conductivity  (S  m- 1 ) 

r  Bruggeman  factor 

Tpc  phase  change  characteristic  time  (s) 
r  shear  stress  tensor  (N  m-2) 

0  electrical  potential  (V) 

6)i  production  rate  of  ith  species  due  to  electrochem¬ 

ical  reactions  (kgm-3  s-1) 

A  Species  concentration  (mol  m- 3 ) 

Subscript 

0  reference  state 

a  anode  side  of  the  membrane 

b  bulk  property 

c  cathode  side  of  the  membrane 

ch  channel 

1  liquid 

m  membrane  property 

p  pore  property 

rxn  reaction 

s  solid 

Superscript 

sat  saturation 


rent  CFCE  tools  is  required  because  these  models  lack  either 
generality  and/or  rational  foundations  [10]. 

In  this  paper  and  its  companion  Part  2  [11]  we  report  on 
the  development  of  a  3D  simulation  code  for  PEMFC  and  its 
validation  using  spatially  resolved  experimental  data  of  local 
current  density  and  water  mass  in  the  MEA,  with  the  pur¬ 
pose  of  establishing  the  bounds  of  validity  of  state-of-the-art 
CFCE  methodology  and  identifying  critical  issues  to  improve 
fidelity  and  reliability.  The  paper  begins  with  a  description  of 
the  governing  equations,  including  the  phenomenological  mem¬ 
brane  model,  and  the  coupling  of  these  equations,  followed  by 
a  description  of  the  computational  domain  and  boundary  condi¬ 
tions.  In  Section  4,  the  experimental  setup  and  post-processing 
of  the  computational  results  are  first  discussed.  Detailed  com¬ 
parison  of  the  numerical  predictions  with  spatially  resolved 
experimental  data,  as  well  as  a  parametric  study  and  analysis 
are  presented  in  Part  2. 

2.  Mathematical  formulation 

2.7.  Governing  equations 

The  CFCE  framework  combines  established  CFD  method¬ 
ology,  which  solves  the  discretized  conservation  equations  for 
mass,  momentum,  species,  and  energy,  with  additional  conser¬ 
vation  equations  needed  to  account  for  electrochemical  reaction 
kinetics  and  the  transport  of  charged  species,  which  in  gen- 
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eral  appear  as  the  potentials  in  the  electrolyte  phase  and  the 
electric  conductor  phase  respectively.  In  addition,  in  a  PEMFC 
two-phase  flow  and  phase  change  of  water  are  likely  to  occur, 
which  necessitates  another  set  of  model  conservation  equations 
to  describe  the  transport  of  liquid  water.  The  multi-phase,  multi¬ 
dimensional  non-isothermal  CFCE  framework  developed  here 
is  based  on  the  following  assumptions: 

(1)  Ideal  gas  for  all  gas  species. 

(2)  No  deformation  in  all  parts  of  the  cell  (no  swelling/shrinking 
of  the  membrane,  no  deformation  of  the  GDE  under  the  land 
area  due  to  compression). 

(3)  Phase  equilibrium  of  water  with  the  electrolyte.  This  allows 
the  use  of  membrane  sorption  isotherm  using  the  water 
activity  at  the  membrane  boundaries. 

(4)  The  unit  cell  operates  at  steady  state. 

Additional  assumptions  required  in  modeling  specific  trans¬ 
port  processes  in  components  will  be  discussed  subsequently.  A 
commercial  CFD  code,  CFD-ACE+,  is  employed  to  solve  the 
complete  set  of  conservation  equations  for  a  unit  cell  geometry. 

2.1.1.  Conservation  of  mass  and  momentum 

The  volume- averaged  mass  and  momentum  conservation 
equations  in  a  porous  media  are: 

fsp)  +  S/{epu)  =0  (1) 

Ot 

and 

3  _  _  £2LLU 

—  (spu)  +  V(spuu)  =  —  eV p  +  V(£t)  +  sB  -| - (2) 

3 1  K 


by 

E  vo  = 0  (6) 

j=l  Arxn 

In  a  PEMFC,  charge  transport  consists  of  protonic  and  electronic 
currents,  and  Eq.  (6)  can  be  written  as: 

-V/h+  =  Vie-  =  Jt  (7) 

The  transfer  current  density  is  non-zero  only  in  the  region  where 
electrochemical  reactions  take  place.  Introducing  an  electrical 
potential  for  each  charged  species  and  relating  the  potential  to 
current  density  using  Ohm’s  law,  we  have 

V(o-mV0m)  =  -V(ctsV0s)  =  jj  (8) 


where  0m  and  0S  are  the  electric  potentials  of  proton  and  elec¬ 
tron,  respectively.  The  transfer  current,  Jj ,  can  be  described 
under  normal  conditions  by  the  Butler- Volmer  equation: 


JTJ 


JOJ 


YlLMk,reffkJ 


exp 


aa  JF 

RUT 


-exp 


«c  ,jF 
- - — ] 

RuT 


N 

k=  1 


(9) 


where  [ A &]  represents  the  average  interfacial  molar  concentra¬ 
tion  of  the  Mi  species  and  the  overpotential  for  each  reaction  is 
defined  as  r)  =  0S  —  0m. 

For  an  electrochemical  reaction  having  the  generalized  form 


NG 

E4A«'  ±  e 


i= 1 


ng 

^  'aijAj  for  j =  1,2,...,  A^steps 
i=  1 


(10) 


2.1.2.  Conservation  of  non- charged  species 

The  mass  conservation  equations  for  the  individual  gas-phase 
species,  i  =  1 ,  . . .,  Ng,  may  be  written  as 

epYi )  +  V(spuYi)  =  VJi  +  cot  (3) 

ot 

where  in  general,  the  species  diffusion  flux  is  given  by  the 
Stefan-Max  well  equation: 

Ji  =  pDjVYj  +  fj-D.VM  -  pY^DjVYj  -  ff-^OjYj 

j  j 

(4) 

The  effective  mass  diffusion  coefficients  of  species  i  within  the 
porous  medium  DL  is  related  to  the  porosity  via  the  Bruggeman 
relation: 

A  =  A,oef  (5) 

where  the  empirical  exponent  §  is  usually  taken  equal  to  1.5. 
The  actual  values  used  in  this  work  are  discussed  in  Section  4. 

2.1.3.  Conservation  of  charged  species  and 
electrochemical  reactions 

Electro-neutrality  dictates  that  the  sum  of  all  current  flows 
equals  to  zero,  and  conservation  of  charged  species  is  thus  given 


the  transfer  current  density  is  related  to  gas  species  concentration 
by  the  so-called  reaction-diffusion  balance  equation: 


v  steps 


7=1 


Yi  -  YPj 


(ID 


where  Ypj  denotes  the  species  mass  fraction  in  the  pore-fluid, 
while  Y(  denotes  the  mass  fraction  at  the  pore-fluid/catalyst  inter¬ 
face.  The  mass  transfer  in  the  pore  level  is  assumed  to  be  limited 
by  the  diffusion  from  the  pore  to  the  catalyst  sites,  hence  an 
average  pore  size  is  used  for  8. 


2.1.4.  Conservation  of  energy 

The  governing  equations  used  in  the  present  work  differ  from 
those  reported  in  Mazumder  and  Cole  [5,6]  mainly  in  the  energy 
equation  and  the  transport  of  liquid  water  (saturation).  In  [5]  the 
energy  equation  is  expressed  and  solved  in  terms  of  enthalpy, 
which  consists  of  enthalpy  of  formation  and  sensible  heat.  In 
the  present  study,  sensible  heat  is  the  primary  variable  for  the 
energy  equation.  While  combining  enthalpy  of  formation  and 
sensible  heat  in  the  energy  equation  is  a  common  practice  when 
dealing  with  reacting  flows  such  as  combustion,  this  formula¬ 
tion  is  not  as  robust  for  an  electrochemical  system  involving 
half-cell  reactions.  The  electrical  power  generated  in  each  half¬ 
cell  reaction  in  this  formulation  should  appear  as  a  sink  term  in 
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the  conservation  equation;  however,  this  is  not  feasible  because 
the  charged  species  (H+,  e-)  and  electrically  neutral  species 
(H2,  O2,  and  H2O)  use  different  reference  points,  which  makes 
evaluation  of  the  Gibbs’  free  energy  for  each  half-cell  reaction 
difficult.  This  ambiguity  can  be  circumvented  by  considering 
only  the  conservation  of  the  sensible  heat  of  the  neutral  species. 
All  potential  losses  then  appear  as  heat  source  terms  in  the  energy 
equation,  and  the  available  electrical  power  is  the  outcome  of 
the  computation.  The  conservation  of  energy  is  written  as: 


a 

—  ((1  £)pshs  +  eph)  +  V  •  (spuh) 
ot 


—  V  •  q  £T  :  V u  £ 


dp 
d  t 


+  h 


i  •  i 


o 


+ 


(12) 


where 


ng 

q  =  XS/T  +  YJJi(hi-gi)  (13) 

i= 1 


The  quantity  —  h(S/  E)eff?7  *s  the  electrical  work 

done  during  an  electrochemical  reaction.  Note  that 


NG  ^ 

Y JiQii  -  gi)  +  Jt 
i—\ 


NG 


eff 


=  Y?iTsi+rt  v  71  (14) 


i=l 


eff 


where  the  first  term  on  the  right  hand  side  represents  reversible 
heat  generation  and  the  second  term  represents  irreversible  heat 
generation.  These  losses  are  in  addition  to  the  Joule  heating  term, 
which  also  causes  irreversible  heat  generation  in  the  cell. 


2.1.5.  Transport  of  liquid  water 

The  transport  of  liquid  water  is  described  in  terms  of 
conservation  of  liquid  saturation,  taking  into  account  evapora¬ 


tion/condensation,  capillary  diffusion  and  gravity: 


d 

—(sdpis)  +  V(sdpus)  =  V(£dpiDcVs) 
ot 

A(1  -  s)k(p\  -  pg)  A  . 

-V  (  - — g  +  m\ 


The  evaporation/condensation  rate  is  expressed  as: 


m\ 


£dP(Yu20  -  ^H2o)/Tpc 

0;  if  s  <0  and  fn\  >  0 


(15) 


(16) 


Eq.  (15)  is  applied  for  flows  in  porous  media  as  well  as  ordinary 
fluid  flow  for  liquid  water  saturation,  even  though  in  the  ordinary 
fluid  phase,  liquid  saturation  loses  its  physical  meaning  as  vol¬ 
ume  fraction  of  liquid  water  in  the  pore.  The  capillary  pressure 
formulation  used  in  [6]  is  adopted  with  a  slight  modification  in 
the  rate  of  phase  change. 


2.2.  Membrane  model 


In  Mazumder  and  Cole  [5]  the  transport  of  non- vapor  water 
is  modeled  as  liquid,  and  water  movement  due  to  the  electro- 
osmotic  drag  (EOD)  is  included  as  a  convective  term  in  the 
conservation  equation  for  liquid  water.  One  drawback  with  this 
approach  is  that  it  fails  to  account  for  water  transport  inside 
the  membrane  phase  when  the  vapor  phase  is  at  undersatu¬ 
rated  conditions.  In  the  present  study,  a  separate  model  for  water 
transport  in  the  membrane  phase  is  adopted.  The  phenomeno¬ 
logical  model  describes  the  water  movement  due  to  the  EOD 
and  diffusion  and  is  solved  only  for  the  membrane  part  of  the 
unit  cell.  The  implementation  of  this  model  within  the  CFD 
framework  requires  special  treatment  for  mass  transfer  on  the 
membrane-catalyst  layer  interface  to  ensure  continuity  of  water 
species. 

The  mechanisms  for  transport  of  water  in  the  electrolyte 
phase  differ  from  those  in  vapor  phase.  The  dominant  water 


Table  1 

Conservation  equations  solved  in  the  comprehensive  numerical  simulation 


Convection 

Diffusion 

Source 

Mass 

V(spu) 

- 

m 

Momentum 

V(spuu) 

—VP  +  V(ef) 

-  (e2%)  n 

v(wr)  +  ^//cp>fvr 

NG 

li-l 

Energy 

sCvpuW(T) 

hn  +  Yjiihj  -  gp 

(|)eff  +  E  +  5h  +  e 

Species 

V(epuYi) 

i 

V7;a 

f 

Potential 

- 

V(cr kV<Pk) 

Saturation 

V(s&pus) 

VKiPi  Dc  V.v] 

p  (  mwpsat  y  \  b 

£&x\mP  /w) 

Water  content 

v(#) 

-$-v(D>yx) 

a  Ji  =  pDiVYi  +  §DiVM 

-pY^DjWYj- 

pYiAM\^n_y_ 

M  / 

j  j 


b  V  =  max([/d,  U),  Ui=  2s.,  L  = 


E- 


C  vol(ic) 


1/3 


Nc 
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transport  mechanisms  in  the  membrane  include  the  EOD  asso¬ 
ciated  with  the  movement  of  charged  species,  diffusion  driven 
by  gradient  of  chemical  potential,  and  hydraulic  permeation  due 
to  pressure  gradients.  In  the  present  work,  the  contribution  due  to 
pressure  gradient  is  assumed  negligible  based  on  experimental 
observation  that  the  membrane  permeability  is  extremely  low  for 
gas  and  liquid.  The  flux  of  water  inside  the  membrane  is  replaced 
by  the  net  water  flux  given  by  a  phenomenological  model  of  a 
form  similar  to  Eq.  (4),  consisting  of  an  EOD  term  and  a  dif¬ 
fusion  term.  The  flux,  expressed  in  terms  of  water  content  X 
(number  of  water  molecules  absorbed  per  sulphonic  acid  site), 
is  thus  given  as: 


h  =  ~i  -  ^DkVX 
F  Mm 


(17) 


Both  the  drag  coefficient  and  water  diffusion  coefficient  are 
functions  of  water  content  and  temperature: 


n d  =  nd(X,  T ), 


(18) 


and 


DX  =  Dx(X,  T)  (19) 

In  the  present  study  the  water  transport  in  the  membrane  is  solved 
separately  from  the  gas  species  equation.  On  the  membrane 
boundaries,  the  flux  of  water  calculated  using  (17)  enters  the  cat¬ 
alyst  layer  domain,  while  the  water  content  boundary  condition 
on  the  membrane  boundaries  is  obtained  through  an  adsorp¬ 
tion  isotherm  for  water  in  the  electrolyte  phase.  The  adsorption 
isotherm  is  a  phase  equilibrium  condition  that  relates  the  water 
content  to  water  activity  and  temperature  in  the  vapor  phase  on 
the  membrane  boundary: 

X  =  X(a,  T)  (20) 


where  the  water  activity  (a)  corresponds  to  the  relative  humidity 
of  water  vapor.  The  ionic  conductivity  of  the  electrolyte  is  gen¬ 
erally  considered  a  function  of  water  content  and  temperature: 

a  =  a(X,  T)  (21) 


Temperature 

(Energy)  Latent  heat 


<J>  source 

Electrical  ^ 
potential  y 

lt*mhnn»*  rr 


Fig.  1.  Coupling  of  conservation  equations  in  PEMFC  simulation. 


heat  in  the  in-plane  and  the  through-plane  directions.  The  cat¬ 
alyst  layer  is  a  composite  material  that  is  constructed  by  two 
networks  (electrolyte  and  carbon  black  respectively),  which 
shows  some  degree  of  anisotropy  due  to  fabrication  processes 
and  also  possibly  due  to  preferred  material  orientation.  The 
bipolar  plate  that  is  made  of  graphite  powder  and  resin  binder 
has  some  flaky  structure  that  makes  its  transport  properties 
anisotropic.  Among  all  these  materials,  the  anisotropy  in  the 
GDL  is  believed  to  be  the  most  influential  because  of  its  fibrous 
structure  as  well  as  the  fact  that  it  is  subject  to  discontinuous 
property  change  on  one  boundary,  i.e.  on  the  side  that  is  in  con¬ 
tact  with  the  bipolar  plate  [17].  Furthermore,  part  of  the  GDL 
on  the  GDL-bipolar  plate  interface  is  under  compression  and  the 
transport  through  such  regions  is  complicated  [18,19].  Owing  to 
the  complexity  of  coupled  transport  phenomena,  analysis  on  the 
effects  due  to  anisotropy  in  transport  properties  of  the  materials 
in  a  PEMFC  is  best  assisted  by  employing  a  simulation  tool. 

The  transport  property  for  a  thermal  or  electrical  conduction 
equation  takes  the  form  of  a  tensor: 


'xx 

kXy 

kxz 

fx 

kyy 

kyz 

fx 

kZy 

kzz 

(22) 


Substituting  for  the  water  flux  in  Eq.  (3)  with  Eq.  (17)  and 
using  transport  properties  of  Eqs.  (18),  (19)  and  (21),  and  the 
phase  equilibrium  relationship  of  (20),  the  transport  of  water 
in  the  membrane  is  solved  and  coupled  iteratively  with  the 
domain  outside  of  the  membrane.  Numerical  analysis  and  behav¬ 
ior  of  this  formulation  is  discussed  in  Sui  and  Djilali  [12]  and 
Mazumder  [13].  The  properties  reported  by  Springer  et  al.  [14], 
cf.  Appendix  B,  are  used  for  the  calculation  unless  otherwise 
noted.  More  details  regarding  the  phenomenological  model  can 
be  found  in  Janssen  [15],  and  a  critical  examination  and  analy¬ 
sis  of  the  transport  models  for  polymer  membranes  is  given  in 
Fimrite  et  al.  [16]. 

2.2.7.  Anisotropic  transport  properties  for  solid  and  porous 
media 

Most  of  the  components  of  a  PEMFC  have  anisotropic  trans¬ 
port  properties.  The  gas  diffusion  layer  (GDL)  in  particular 
shows  significant  differences  in  the  transport  of  electricity  and 


For  a  conventional  plate-and-frame  type  of  membrane  elec¬ 
trode  assembly  (ME A)  with  coordinates  X  and  Z  as  the  in-plane 
direction  and  Y  as  the  through-plane  direction,  we  assume  the 
transport  is  limited  in  both  planes,  i.e.  only  diagonal  terms  in  the 
matrix  remain.  For  a  PEMFC  problem,  the  anisotropy  should 
appear  in  the  conservation  of  species  (diffusivity),  momentum 
(permeability),  energy  (heat  conductivity)  and  potential  (elec¬ 
trical  conductivity). 

2.3.  Coupling  of  transport  equations 

Table  1  summarizes  the  governing  equations  implemented 
and  solved  in  the  CFD  code.  One  of  the  most  challenging  aspects 
from  a  computational  view  point  is  the  coupling  between  the 
various  transport  phenomena  within  a  PEMFC  as  illustrated  in 
Fig.  1 .  Each  circle  represents  a  conservation  equation  and  the 
outward-pointing  arrows  represent  the  effects  of  the  conserva¬ 
tion  equation  in  question  on  other  equations.  The  conservation 
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equation  of  gas  species,  water  in  liquid  and  solid  (membrane) 
phases,  and  energy  are  closely  coupled  with  all  other  conserva¬ 
tion  equations.  Among  all  the  processes  considered,  the  transport 
of  water,  which  exists  in  the  system  in  the  forms  of  vapor,  liquid, 
and  absorbed  in  solid,  is  central  to  all  other  coupled  transport 
processes.  Unsaturated  water  in  the  vapor  phase  affects  local  rel¬ 
ative  humidity  near  the  membrane,  especially  on  the  anode  side, 
and  thus  affects  water  transport  across  the  membrane  as  well 
as  electrical  properties  of  the  membrane,  which  in  turn  impacts 
the  solution  for  the  electric  potential.  The  transport  of  liquid 
water  in  the  porous  media  affects  mass  transport  of  gas  species, 
while  liquid  water  in  the  gas  channels  changes  the  pressure  field 
and  may  alter  flow  distribution  in  a  unit  cell.  Transport  of  water 
across  the  membrane  is  a  key  phenomenon  that  links  the  trans¬ 
port  processes  between  the  anode  and  cathode  sides.  Therefore 
modeling  of  water  transport  in  the  membrane  phase  is  crucial  to 
the  simulation  capability  development  for  PEMFC. 

The  electrical  potential  plays  two  roles  in  the  coupling.  On 
the  one  hand  it  can  be  considered  as  representing  conservation 
of  charged  species  (protons  and  electrons),  and  on  the  other 
hand  it  can  be  viewed  as  part  of  the  conservation  of  energy,  in 
which  the  electrical  power  converted  from  the  electrochemical 
reaction  manifests  in  a  potential.  The  electrical  potentials  are 
conveniently  solved  separately  from  the  species  equations  to 
avoid  the  implicitness  of  the  potential  difference  that  appears  in 
the  Butler- Volmer  equation  used  for  the  reaction  rate  in  the  gas 
species  equations.  This  effectively  decouples  the  two  half-cell 
reactions  in  the  anode  and  cathode,  and  as  a  result  convergence 
in  the  gas  species  is  not  as  robust  as  that  in  the  conventional 
treatment  of  chemical  reactions  in  CFD.  While  it  is  possible  to 
include  the  potential  into  the  energy  equation,  the  difficulty  is 
that  the  chemical  energy  and  the  electrical  energy  have  different 
reference  points. 


3.  Numerical  implementation 

To  ensure  continuity  of  water  transfer  across  the  membrane- 
catalyst  layer  interface,  the  water  content  value  on  this  interface 
is  calculated  by 


U  Yw  n^(Xw)iw  Xw  Xm 

-pD - - -  =  - - - PmDkfrc) - * -  (23) 

8C  F  8m 

Fig.  2  shows  the  relation  of  the  variables  in  Eq.  (17)  in  the 
vicinity  of  the  membrane  boundary.  The  proper  implementation 
of  the  membrane  model  into  CFD-ACE  was  checked  against  a 
ID  solver  for  Eq.  (1),  cf.  Sui  and  Djilali  [12],  with  prescribed 
current  density  and  boundary  condition  values  on  the  membrane, 
see  Fig.  3.  The  predicted  polarization  curve  obtained  using  the 
code  with  the  membrane  water  transport  model  as  well  as  the 
curve  predicted  with  the  code  described  in  Mazumder  and  Cole 
[5]  are  compared  to  the  experimental  data  of  Ticianelli  et  al. 
[20]  in  Fig.  4,  showing  a  clear  improvement  in  the  predictions 
with  the  implementation  of  the  membrane  water  transport  model. 
These  predictions  were  made  by  setting  all  model  parameters  to 
correspond  to  those  given  in  [20]. 


Membrane 

X 


Fig.  2.  Boundary  condition  on  the  membrane/catalyst  layer  interface. 


Fig.  3.  Numerical  validation  of  the  membrane  model  implemented  in  the  sim¬ 
ulation. 


Polarization  curve  for  323 K  and  1/1  atm 


Current  density  (A/cm2) 

Fig.  4.  Comparison  of  polarization  curves  by  experiment  and  2D  predictions 
with  and  w/o  Springer  model  for  membrane.  Springer  model  apparently  predicts 
the  current  more  close  to  measurement  because  the  water  transport  across  the 
membrane  is  better  accounted  for. 


Table  2 


Summary  of  properties  and  parameters  used  for  the  baseline  calculation 


Unit 

Bipolar 

plate 

Anode  gas 
channel 

Anode  GDL 

Anode  catalyst  layer 

Membrane 

Cathode  Catalyst  layer 

Cathode  GDL 

Cathode  gas 
channel 

Coolant 

channel 

Dimension 

m 

- 

2.5  x  10-4a 

2  x  lO"4 

3.5  x  10“5 

5.4  x  10-5 

3.5  x  10“5 

2.4  x  10"4 

3.9  x  10~5a 

3.5  x  10“5a 

Porosity 

- 

- 

0.6 

0.6 

0.1 

0.6 

0.6 

- 

- 

Average  pore 

m 

- 

- 

2.5  x  10"6 

1  x  10"7 

1  x  10“9 

1  x  10“7 

2.5  x  10"6 

- 

- 

size 

Bruggeman  £ 
Permeability 

m2 

- 

- 

1 

1  X  10“13 

2 

1  x  10“13 

13 

1  x  10“21 

2 

1  x  10“13 

1 

1  X  10“13 

- 

- 

Thermal 

W  m-1  K-1 

- 

20 

20 

20 

20 

20 

- 

- 

conductivity 

Electrical 

(£2  m)  1 

- 

200 

80 

- 

80 

200 

- 

- 

conductivity 

Protonic 

(Qm)-1 

- 

- 

- 

5 

Springer 

5 

- 

- 

- 

conductivity 

Diffusion 

m2  s  1 

- 

Sc  =  0.7 

Sc  =  0.1 

& 

ii 

o 

Ci 

Springer 

& 

ii 

o 

C 

5c  =  0.7 

ii 

o 

Ci 

- 

coefficient 

Density 

kgm  3 

1600 

IG 

IG 

IG 

1980 

IG 

IG 

IG 

Viscosity 

kgm-1  s-1 

- 

MKT 

MKT 

MKT 

- 

MKT 

MKT 

MKT 

Reaction 

- 

- 

- 

HOR 

- 

ORR 

- 

- 

- 

Coefficients  for 

reaction 

70  =  Am-3, 
S!V=m-' 

aa  =  0.5,  ac  =  0.5, 
jo  =  1  x  10-9,  yh2  =  1, 
S/V=  1000 

aa  =  1,  ac  =  1.2, 

jo  =  5  x  10“6,  yo2  =  1, 

S/V  =  1 000 

IG,  ideal  gas  [21];  MKT,  mixed  kinetic  theory  [21];  HOR,  hydrogen  oxidation  reaction;  ORR,  oxygen  reduction  reaction;  Sc,  Schmidt  number.  Unless  otherwise  stated,  all  other  data  obtained  from  [25]. 
a  Hydraulic  diameter. 
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Fig.  5.  Computational  results  of  the  baseline  case. 


3.1.  Computational  domain  and  boundary  conditions 

The  generic  computational  model  was  applied  to  simulate  a 
straight  channel  unit  cell  configuration  related  but  not  identical 
to  the  Ballard  Mk902  hardware,  cf.  Fig.  5.  Properly  designed 
and  optimized,  such  straight  channel  configurations  can  offer 
a  good  balance  between  various  thermo-fluid  and  manufactur¬ 
ing  requirements,  and  are  becoming  more  common  in  industrial 
stacks.  The  computational  model  corresponds  to  the  experimen¬ 
tal  cell  which  consists  of  an  array  of  evenly  spaced  straight 
channels  with  a  length  of  0.6  m.  The  computational  domain 
includes  the  segment  bounded  by  the  center  lines  of  two  con¬ 
secutive  land  areas  with  the  full  length  of  the  unit  cell,  which 
includes  the  ME  A,  gas  and  coolant  channels,  all  surrounded 
by  the  bipolar  material.  All  cases  reported  are  with  anode  and 
cathode  gases  flowing  in  opposite  directions  (‘counter-flow’) 
and  coolant  flowing  in  the  same  direction  as  the  cathode  flow. 
An  adiabatic  boundary  condition  is  applied  to  all  surfaces  of 
the  computational  domain  except  for  the  openings  of  the  flow 
channels,  which  have  prescribed  temperature  conditions. 

3.2.  Post-processing  of  3D  data 

The  experimental  data  provides  spatial  resolution  along  the 
channel  direction  but  the  distributions  are  essentially  averaged 
in  the  span  wise  direction.  In  order  to  allow  direct  comparison, 
the  3D  computational  results  are  integrated  in  the  appropriate 
dimensions.  For  the  channel  flows,  the  flow  quantities  are  inte¬ 
grated  at  each  location  over  the  channel  cross-section  using  the 
following  equations: 

Mass  flow  rate: 

rhi(z)  =  J  puYidAch ,  (24) 

Liquid  water  in  gas  channel: 

mw(z)  =  /  pwusdAch,  (25) 


Bulk  fluid  temperature: 

f  puCpT  dAch 
b"  /  puCp  dAch 


(26) 


The  water  balance  is  calculated  using  the  inlet  and  outlet  values 
of  water  mass  flow  rates: 


_  (^H20,QUT  -  WH20,IN)ai10c|e 
m  Hi  O,  generated 

_  (WHrO.OUT  ~  WH2Q,IN)anocie 
(wo2JN  —  rn  o2 , 0  UT  )2Mh2q/M  02 


(27) 


The  average  current  density  is  calculated  on  the  outer  surface  of 
the  cathode  current  collector: 


m  = 


f  idAc 


(28) 


Membrane  resistance: 


R(z)  = 


JjQ?) 


dAc 


(29) 


The  water  weight  per  puck  is  calculated  by  adding  the  water  in 
the  electrolyte  phase  and  liquid  water  in  the  pore: 


^H20,y(ch°p)  —  J  Pinion omerk  d Vj  +  r jSpn2o, 

j  =  GDL,  CL,  PEM 


(30) 


4.  Results  and  discussion 

4.1.  Parameter  and  properties  used  for  the  baseline  case 

Table  2  lists  the  parameter  values  and  properties  used  for 
the  baseline  calculations.  The  calculation  method  of  standard 
transport  properties  (density,  diffusion  coefficient,  viscosity, 
thermal  conductivity)  for  gas  species  are  given  in  [21].  The  per¬ 
meability  and  porosity  and  the  Bruggeman  coefficient  for  the 
membrane  are  chosen  to  effectively  preclude  permeation  across 
the  membrane  due  to  a  pressure  difference  across  the  mem¬ 
brane.  The  electrical  conductivity  of  the  GDL  is  isotropic  for 
the  baseline  case  but  calculations  using  anisotropic  conductiv¬ 
ities  are  possible.  The  protonic  conductivity  of  the  electrolyte 
phase  in  the  catalyst  layers  is  set  to  a  constant  for  the  base¬ 
line  for  simplicity.  The  transfer  coefficient  for  oxygen  reduction 
reaction  (ORR),  ac ,  is  the  only  control  parameter  anchored  to 
fit  the  polarization  data.  Once  anchored  at  a  particular  con¬ 
dition,  the  same  value  is  used  to  generate  the  performance 
curve. 

4.2.  Baseline  conditions 

Fig.  5  shows  some  computational  results  of  the  baseline  case. 
In  the  3D  model,  one  can  see  gradient  in  species  mass  fraction 
and  current  flow  in  the  direction  between  the  channel  area  and  the 
land  area.  Although  these  gradients  may  have  significant  impact 
on  cell  performance  and  durability,  such  gradients  are  not  the 
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Fig.  6.  Mass  flow  rate  of  gas  and  liquid  in  the  gas  channels. 


Fig.  8.  Current  density  and  membrane  resistance. 


focus  of  the  present  study;  readers  are  referred  to  the  recent 
work  of  Freunberger  et  al.  [22]  for  a  discussion  of  this  issue.  In 
the  present  study,  we  focus  on  validation  based  on  the  gradients 
in  the  axial  direction  (along  the  channel)  to  coincide  with  the 
available  experimental  data.  Figs.  6-9  show  the  computational 
results  of  the  baseline  case  plotted  in  the  axial  direction  of  the 
unit  cell  with  the  cathode  inlet  (anode  outlet)  located  Z-  0  and 
the  cathode  outlet  (anode  inlet)  at  Z  =  1 .  The  transport  proper¬ 
ties  and  dimensions  used  in  the  computation  are  summarized  in 
Table  2.  Operating  conditions  of  the  unit  cell  for  the  baseline 
case  are:  7=1  A  cm-2;  F  =  3bars,  RH  =  55%,  inlet  tempera¬ 
ture  =65.1  °C,  and  stoichiometric  ratio  =  1.71  for  the  cathode; 
andP  =  3.2  bars,  RH  =  19%,  inlet  temperature  =  73.9  °C,  and  sto¬ 
ichiometric  ratio  =  1.54  for  the  anode.  Fig.  6  shows  the  mass 
flow  rates  of  gas  species  and  liquid  water  along  the  channel 
calculated  using  Eqs.  (24)  and  (25),  respectively.  It  should  be 
noted  that  the  mass  flow  rates  in  the  anode  channel  (right  verti¬ 
cal  axis  in  Fig.  6)  are  roughly  an  order  of  magnitude  less  than 
in  the  cathode  side  (left  axis  in  Fig.  6).  The  fuel  and  oxidant 


flow  rates  appear  to  decrease  almost  linearly  along  the  chan¬ 
nel,  which  indicates  the  local  flow  rate  into  the  GDL/channel 
interface  is  nearly  constant  at  any  location  of  the  channel.  How¬ 
ever  this  flow  redistributes  inside  the  MEA  by  diffusion  based 
on  local  potential  distribution  under  the  land  area.  The  local 
consumption  rates  of  the  oxidant  and  the  fuel  thus  become  less 
closely  related  to  local  current  density  distribution.  The  com¬ 
mon  assumption  made  for  along-the-channel  type  of  models, 
see  e.g.  [23,24],  whereby  local  reactant  consumption  is  related 
to  local  current  density  is  therefore  not  valid.  The  water  vapor  in 
the  anode  channel  shows  an  increase  near  the  inlet  portion  fol¬ 
lowed  by  a  decrease  downstream.  The  increase  of  water  vapor 
flow  close  to  the  inlet  is  mainly  due  to  back  diffusion  from  the 
cathode  outlet,  which  is  likely  fully  humidified.  As  the  anode 
channel  is  gradually  humidified,  the  EOD  begins  to  take  effect 
and  drags  increasingly  more  water  into  the  cathode  side;  thus  a 
pattern  of  net  water  recirculation  forms  within  the  unit  cell.  The 
water  vapor  flow  rate  in  the  cathode  first  increases  slightly  near 
the  inlet,  and  then  stabilizes  as  the  cathode  stream  becomes  fully 
saturated.  Liquid  water  only  appears  in  the  cathode  channel  for 


Fig.  7.  Relative  humidity  in  the  gas  channels. 


Fig.  9.  Temperature  profiles  in  the  gas  channels. 
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the  baseline  case.  As  the  cathode  flow  proceeds  downstream,  it 
takes  up  more  product  water  from  the  ORR,  but  near  the  cathode 
outlet  a  slight  decrease  in  the  liquid  flow  rate  is  observed  due 
to  back  diffusion  of  water  through  the  membrane  to  the  anode. 
Fig.  7  shows  the  RH  calculated  based  on  bulk  fluid  tempera¬ 
ture  and  vapor  flow  rate  in  the  anode  and  cathode  channels.  The 
water  diffusion  across  the  membrane  is  primarily  dependent  on 
the  water  activities  on  both  sides  of  the  membrane,  which  is 
closely  related  to  the  RH  in  the  channels.  When  RHA>RHC 
(approximately  Z  <  0. 1  in  Fig.  7),  diffusion  of  water  in  the  mem¬ 
brane  is  in  the  same  direction  of  the  EOD,  while  in  the  rest  of 
the  unit  cell,  water  diffusion  opposes  water  transport  via  the 
EOD.  Near  the  anode  inlet,  where  the  difference  in  RH  between 
anode  and  cathode  is  largest,  back  diffusion  from  the  cathode 
side  dominates  water  transport  across  the  membrane  and  the  net 
transport  of  water  is  from  cathode  to  anode. 

Fig.  8  shows  the  current  density  profile  and  membrane  resis¬ 
tance  along  the  flow  channels.  Two  higher  resistance  regions 
can  be  observed  near  the  anode  and  cathode  inlets  respectively 
as  a  result  of  the  low  RH  values  in  these  regions.  The  higher 
resistance  causes  a  small  current  fall  off  near  the  cathode  inlet 
and  a  larger  drop  near  the  anode  inlet,  where  the  high  resis¬ 
tance  is  compounded  by  lower  oxygen  concentration  on  the 
cathode  side.  Fig.  9  shows  the  bulk  fluid  temperature  profiles 
for  the  reactant  and  coolant  channels.  One  can  see  that  the 
coolant  temperature  varies  almost  linearly  along  the  coolant 
channel,  with  a  temperature  increase  of  about  8  °C.  It  should 
be  noted  that  in  general  serpentine  flow  channel  cells  would  be 
expected  to  exhibit  lower  temperature  differences  due  the  lower 
geometric  aspect  ratio  and  enhanced  heat  transfer  due  to  the 
inherent  cross  and  counter  flow  feature  of  serpentine  arrange¬ 
ments. 

4.3.  Similarity  of  computational  results 

The  straight  channel  unit  cell  geometry  considered  here 
enables  not  only  less  ambiguous  profile  measurement  for  com¬ 
parison  with  numerical  results,  but  is  also  convenient  for 
assessing  the  validity  of  simplified  models  with  reduced  dimen¬ 
sion.  Since  the  length  scale  of  the  cell  along  the  flow  is  much 
larger  than  all  other  dimensions,  the  gradients  in  the  axial  direc¬ 
tion  are  expected  to  be  relatively  smaller.  This  suggests  that 
locally  the  solutions  of  all  variables  are  dictated  by  the  2D  cross- 
section  perpendicular  to  the  axial  direction.  Fig.  10  shows  the 
profiles  of  oxygen  and  water  mass  fraction  as  well  as  current 
density  on  the  GDL/catalyst  layer  interface  at  ten  evenly  spaced 
axial  locations.  For  each  variable  shown  in  Fig.  10,  these  profiles 
can  be  collapsed  into  one  profile  by  normalizing  the  local  pro¬ 
file  with  the  difference  between  the  maximum  and  the  minimum 
value  of  the  profile,  cf.  Fig.  11. 

The  profile  similarity  suggests  that  a  2D  computational 
domain  may  be  suitable  for  obtaining  a  base  solution  for  the  cou¬ 
pled  transport,  and  with  appropriate  scaling  of  the  base  solution, 
approximate  solutions  can  be  obtained  for  a  3D  geometry.  This 
approach  would  considerably  reduce  computational  resource 
requirements  and  is  attractive  in  the  context  of  systematic  para¬ 
metric  studies  and  optimization. 


Fig.  10.  Profiles  at  10  axial  locations  along  the  channel:  (a)  O2,  (b)  H2O  and 
(c)  and  current  density.  The  numbers  in  each  figure  indicate  the  location  from 
cathode  inlet. 
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presented  together  with  sensitivity  and  parametric  analyses  to 
fully  assess  the  reliability  of  the  computational  tool. 
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Appendix  A.  Energy  conservation  in  a  PEMFC  unit  cell 

The  energy  equation  can  be  written  as: 


Fig.  11.  Normalized  profiles  for  O2,  H2O  mass  fraction  and  current  density  on 
GDL/catalyst  layer  interface. 
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Fig.  Al.  Illustration  of  energy  balance  in  a  fuel  cell. 
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V(epuh)  =  V( -X.VT)  -  V 


(Al) 


vi=  1 


where  the  Gibbs  free  energy  can  be  expressed  for  each  species 
(charged  and  neutral)  in  terms  of  enthalpy  and  product  of  tem¬ 
perature  and  entropy  to  become: 


W(spuh)=W(-XWT)-W 


J2  Ot(ht-Tsi))+  J2  CU4>i) 

V  neutral  charged 


(A2) 


The  enthalpy  in  (A2)  consists  of  the  enthalpy  of  formation  and 
sensible  heat: 

T 

hi  =  Ah°fi  +  JcpdT  =  Ahfi  +  ht  (A3) 

To 

Substituting  (A3)  into  (A2)  we  have 


5.  Conclusions 

We  have  presented  the  modelling  framework  and  develop¬ 
ment  of  a  comprehensive  3D  CFD  based  simulation  code  for 
PEMFC,  with  a  focus  on  implementation  and  base  case  sim¬ 
ulations  to  illustrate  the  physical  features,  transport  of  species 
along  the  channel  and  coupling  between  heat  and  mass  transfer 
processes.  Analysis  of  the  base  case  results  reveals  similarities 
in  the  computed  cross-sectional  profiles  along  the  axial  direction 
of  the  channels,  indicating  that  dominant  coupling  of  transport 
phenomena  is  in  the  2D  cross-section,  This  suggest  that  for  unit 
cells  having  parallel  flow  channels,  reduced  dimensional  model 
based  on  appropriate  similarity  variables  might  provide  a  useful 
and  fast  turnaround  alternative  to  full  CFD  for  initial  design  and 
or  optimization  cycles. 

The  implementation  of  a  phenomological  membrane  water 
transport  was  shown  to  yield  improved  predictions  of  the  global 
polarization  curve.  In  Part  2  of  this  work  [1 1]  a  systematic  vali¬ 
dation  of  the  CFD  code  against  spatially  resolved  experimental 
data,  including  water  mass  balances  and  current  distribution,  is 


=  VC-A.VD  -  V  YfJiiA h°u  +  hi  -  Tsi )  +  ^  (7,0,) 

\i=  1  charged  f 

(A4) 

The  species  equation  from  (A4)  in  Table  1  can  be  written  as 

V(epVYi)  +  V/f  =  rhi  (A5) 

Multiplying  the  enthalpy  of  formation  of  each  species  to  (A5), 
we  have 

Ahf'ViepVYj)  +  Ahfs/Oi)  =  rhiAh°{i  (A6) 

Summing  (A6)  for  all  gas  species,  we  have 

V(epi>  £  A h°fJ)  +  V(7,  ^  A h°fJ)  =  J2  ™^hli  (A7> 

neutral 
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Subtracting  (A7)  from  (A4),  we  have 


spuJ2(hi)  =  V(-AVD  -  V  -  Tsj) 


W=1 


+  (A8) 

charged  J 


Expanding  the  electrical  potential  term  in  (A8),  we  have 


v[  z;  (*.•*))  =o'V0+0vo H+ 

y  charged  y 

+(I  V0  +  0V/)e-  = 

With  V/  =  jt,  we  have 

eCpp3V(T)  +  V(AVD  +  5y(Cp,,vr 

i 

N 

jrn  +  5>*a*) 

i=i 


CT 


+  (V/)i7  (A9) 


+ 


. 

V/eff 


(A10) 


Appendix  B.  Properties  of  the 
Springer-Zawodzinski-Gottesfeld  membrane  model 


Water  content  in  the  electrolyte  phase  is  related  to  water 
activity: 

X  =  0.043  +  17.81a  —  39.85a2  +  36a3  (Bl) 


For  the  vapor  phase  on  the  membrane  surface  the  water  activity  is 
equal  to  the  relative  humidity.  The  drag  coefficient  is  expressed 
as  a  linear  function  of  water  content: 

2.5 

nd  =  —k  (B2) 

The  diffusion  coefficient  given  by  Springer  et  al.  is  written  as 


Dw 


k  da  , 
(1  T  sk )  a  dcw 


(B3) 


where  s  =  0.126  is  the  swelling  factor  and  D'  is  fitted  piecewise 
as 


£>'  =  0.257.,  0  <  7.  <  2 

D'  =  0.5  +  0.8125(A  -  2),  2  <  k  <  6  (B4) 

D'  =  3.75  +  0.267(7.  -  6),  6  <  k 


The  protonic  conductivity  of  the  electrolyte  is  given  by 


cr=(0.005139A— 0.00326)exp 


(B5) 


Because  liquid  water  is  treated  as  independently  as  a  different 
“ species ”,  the  energy  due  to  phase  change  appears  as  a  source 
term  in  (A  10).  Including  the  source  term  for  phase  change  and 
adding  any  possible  external  heat  sources,  we  have  the  conser¬ 
vation  of  energy  in  terms  of  sensible  heat  as: 


Membrane  density  pm  =  2000kgm  3  and  equivalent  weight 
Mw  =  1.1  kg  mol-1  are  used  in  the  calculation. 
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The  electrical  energy  that  comes  out  of  the  system  can  be  cal¬ 
culated  using  the  first  law  of  thermodynamics: 


out 


out 
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Fig.  Al  illustrates  the  relationship  among  the  inlet/outlet  of  the 
system  and  heat/work  in  a  PEMFC  system. 
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